library(gtsummary)
library(dplyr)
library(tidyverse)
library(survival)
library(survminer)
library(lubridate)
library(readxl)
library(janitor)
library(ggplot2)
library(tidycmprsk)
library(ggsurvfit)
library(tidymodels)
library(sjPlot)
library(patchwork)
library(bstfun)
library(ComplexUpset)
library(writexl)
library(rms)
library(timeROC)
library(rms)
library(pec)
library(prodlim)
library(riskRegression)
library(plotly)
library(ggtext)
library(WeightIt)
library(cobalt)
library(survey)
db_Kansas = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 - due 5-1-25 (de-identified).xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
filter(!is.na(age_at_bs_ab_initiation)) %>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
alc_day_8 = as.double(alc_day_8),
alc_day_15 = as.double(alc_day_15),
alc_day_22 = as.double(alc_day_22),
alc_c2d1 = as.double(alc_c2d1),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
)
# Write to Excel
#write_xlsx(db_Kansas, "Kansas_elra_data.xlsx")
db_CC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/CC bsAb Consortium Submission CCF 6-11-25.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
alc_day_8 = as.double(alc_day_8),
alc_day_15 = as.double(alc_day_15),
alc_day_22 = as.double(alc_day_22),
alc_c2d1 = as.double(alc_c2d1),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
)
# Write to Excel
#write_xlsx(db_CC, "CC_elra_data.xlsx")
db_HUMC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/BsAb Consortium Data HUMC.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
alc_day_8 = as.double(alc_day_8),
alc_day_15 = as.double(alc_day_15),
alc_day_22 = as.double(alc_day_22),
alc_c2d1 = as.double(alc_c2d1),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
)
# Write to Excel
#write_xlsx(db_HUMC, "HUMC_elra_data.xlsx")
db_MCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data MCC-deidentified.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
)
db_Huntsman = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/Utah_Huntsman_bsAb Consortium Data FINAL_5.4.25.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc)/1000,
baseline_alc = as.double(baseline_alc)/1000,
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
)
# Write to Excel
#write_xlsx(db_Huntsman, "Huntsman_elra_data.xlsx")
db_FHCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/FHCC Bispecifics (5.5.2025).xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
date_of_infection_80 = ymd(date_of_infection_80),
baseline_anc = as.double(baseline_anc),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
)
db_UTSW = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 1.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3")) %>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles)
)
# Write to Excel
#write_xlsx(db_UTSW, "UTSW_elra_data.xlsx")
db_MDACC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/MDACC BITE ASH 2025.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
)
# Write to Excel
#write_xlsx(db_MDACC, "MDACC_elra_data.xlsx")
db_Roswell = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 5.8.2025_NN.xlsx", sheet = "bsAb_DATA")%>%
filter(`bsAb name` %in% c("1","3"))%>%
clean_names()%>%
mutate(
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
baseline_anc = as.double(baseline_anc),
baseline_alc = as.double(baseline_alc),
anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
)
# Write to Excel
#write_xlsx(db_Roswell, "Roswell_elra_data.xlsx")
# Combine all into a single dataframe
combined_db <- bind_rows(
list(
HUMC = db_HUMC,
MCC = db_MCC,
Huntsman = db_Huntsman,
FHCC = db_FHCC,
UTSW = db_UTSW,
MDACC = db_MDACC,
Roswell = db_Roswell,
CC = db_CC,
Kansas = db_Kansas
),
.id = "site"
) %>%
mutate(
prior_bcma_directed_therapy_yes_no = ifelse(!is.na(most_recent_bcma_agent_name),1,0),
bs_ab_name = case_when(
bs_ab_name == 1 ~ "tec",
bs_ab_name == 3 ~ "elra",
.default = NULL
),
max_crs_grade_0_5 = ifelse(is.na(max_crs_grade_0_5), 0,max_crs_grade_0_5), ## Set NAs to 0
max_icans_grade = ifelse(is.na(max_icans_grade), 0, max_icans_grade ) ## Set NAs to 0
)
# Convert all datetime columns to Date
combined_db_clean <- combined_db %>%
mutate(across(where(lubridate::is.POSIXt), as.Date))
# Write to Excel
write_xlsx(combined_db_clean, "elra_tec_data.xlsx")
# Write to Excel
# write_xlsx(combined_db_clean %>% filter(prior_bcma_directed_therapy_yes_no != prior_bcma_directed_therapy), "elra_data_discordant.xlsx")
## Several institutions do not have any cases of elranatamab
# db_UABMC = read_excel(path = "Data/bsAb Consortium Data 5.8.xlsx", sheet = "bsAb_DATA")%>%
# filter(`bsAb name` == "3")%>%
# clean_names()
# db_Stanford = read_excel(path = "Data/De-identified Stanford bsAb Consortium Data Template 4.30.2025.xlsx", sheet = "bsAb_DATA")%>%
# filter(`bsAb name` == "3")%>%
# clean_names()
# db_LCI = read_excel(path = "Data/bsAb Consortium Data LCI.xlsx", sheet = "bsAb_DATA")%>%
# filter(`bsAb name` == "3")%>%
# clean_names()
# db_COH = read_excel(path = "Data/bsAb Consortium Data Template 2.11.2025_COHv2.xlsx", sheet = "bsAb_DATA")%>%
# filter(`bsAb name` == "3")%>%
# clean_names()
# db_MUSC = read_excel(path = "Data/USMM bsAb Consortium Data Template 2.11.2025.xlsx", sheet = "bsAb_DATA")%>%
# filter(`bsAb name` == "3")%>%
# clean_names()
# db_Mayo = read_excel(path = "Data/Bispecific_consortium_20250521 NoPHI.xlsx", sheet = "bsAb_DATA") %>%
# filter(`bsAb name` == 3)%>%
# clean_names()
combined_db <- combined_db %>%
mutate(
# baseline_anc = as.numeric(baseline_anc),
# baseline_hgb,
# baseline_pl_ts,
# baseline_crp_prior_to_bs_ab_initiation,
# baseline_ferritin_prior_to_bs_ab_initiation,
HT_Score =
(baseline_anc <= 1.2) +
(baseline_hgb <= 9.0) +
(baseline_pl_ts >= 76 & baseline_pl_ts <= 175) +
(baseline_crp_prior_to_bs_ab_initiation >= 3.0) +
(baseline_ferritin_prior_to_bs_ab_initiation >= 650 & baseline_ferritin_prior_to_bs_ab_initiation <= 2000) +
(baseline_pl_ts <= 75) + ## Add an extra point for Plt <=75
(baseline_ferritin_prior_to_bs_ab_initiation > 2000), ## Add an extra point for ferritin >2000
# Classify the risk based on the score
HT_Risk = factor(
case_when(
HT_Score >= 2 ~ "HThigh",
baseline_pl_ts <=75 ~ "HThigh", ## Equivalent to at least 2 points
baseline_ferritin_prior_to_bs_ab_initiation > 2000 ~ "HThigh", ## Equivalent to at least 2 points
(baseline_anc <= 1.2) + (baseline_hgb <= 9.0) + (baseline_pl_ts >= 76 & baseline_pl_ts <= 175) >=2 ~ "HThigh",
HT_Score <2 ~ "HTlow",
.default = "NA"
),
levels = c("HTlow", "HThigh")
)
)
write_xlsx(combined_db %>% select(site, patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate, baseline_anc, baseline_hgb, baseline_pl_ts, baseline_crp_prior_to_bs_ab_initiation, baseline_ferritin_prior_to_bs_ab_initiation, HT_Score, HT_Risk), "HT_score.xlsx")
combined_db <- combined_db %>%
mutate(
severe_infection = ifelse(infection_severity_73 %in% 1 | infection_severity_76 %in% 1 | infection_severity_79 %in% 1| infection_severity_82 %in% 1,1,0),
first_severe_infection_date = case_when(
infection_severity_73 %in% 1 ~ date_of_infection_74,
infection_severity_76 %in% 1 ~ date_of_infection_77,
infection_severity_79 %in% 1 ~ date_of_infection_80,
infection_severity_82 %in% 1 ~ date_of_infection_83,
.default = NULL
),
first_severe_infection_type = case_when(
infection_severity_73 %in% 1 ~ number_infection_1_bacterial_viral_fungal,
infection_severity_76 %in% 1 ~ number_infection_2_bacterial_viral_fungal,
infection_severity_79 %in% 1 ~ number_infection_3_bacterial_viral_fungal,
infection_severity_82 %in% 1 ~ number_infection_4_bacterial_viral_fungal,
.default = NULL
)
)
theme_gtsummary_compact()
variables <- c("elra", "Age", "Female sex", "ECOG 2+", "EMD", "Hgb", "LDH", "Prior BCMA")
data1 <- combined_db %>%
drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
transmute(
Age = age_at_bs_ab_initiation,
"Female sex" = sex_0_male_1_female,
elra = ifelse(bs_ab_name == "elra", 1, 0),
"ECOG 2+" = ifelse(ecog_ps_at_bs_ab_initiation_0_5 >=2,1,0),
EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
Hgb = as.double(baseline_hgb),
LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
"Prior BCMA" = prior_bcma_directed_therapy_yes_no
)
uv_tab <- tbl_uvregression(
data1[c(variables)],
method = glm,
y = "elra",
method.args = list(family = binomial),
exponentiate = TRUE
) %>%
sort_p()
mv_tab<-glm(elra ~ Age + `Female sex` + `ECOG 2+` + EMD + Hgb + LDH + `Prior BCMA`, data=data1, family=binomial) %>%
tbl_regression(exponentiate=TRUE)
tbl_merge(list(uv_tab, mv_tab), tab_spanner = c("**Univariate**", "**Multivariable**"))
| Characteristic |
Univariate
|
Multivariable
|
|||||
|---|---|---|---|---|---|---|---|
| N | OR1 | 95% CI1 | p-value | OR1 | 95% CI1 | p-value | |
| Hgb | 371 | 1.12 | 1.01, 1.25 | 0.029 | 1.14 | 1.02, 1.28 | 0.026 |
| EMD | 371 | 1.69 | 0.95, 2.96 | 0.069 | 1.71 | 0.96, 3.04 | 0.068 |
| Age | 371 | 1.02 | 1.00, 1.04 | 0.087 | 1.01 | 0.99, 1.04 | 0.3 |
| Female sex | 371 | 1.43 | 0.93, 2.21 | 0.11 | 1.43 | 0.92, 2.24 | 0.11 |
| Prior BCMA | 371 | 0.76 | 0.49, 1.17 | 0.2 | 0.78 | 0.50, 1.23 | 0.3 |
| ECOG 2+ | 371 | 1.19 | 0.75, 1.88 | 0.4 | 1.27 | 0.78, 2.07 | 0.3 |
| LDH | 371 | 1.00 | 1.00, 1.00 | 0.8 | 1.00 | 1.00, 1.00 | >0.9 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||||||
IPTW.data <- combined_db %>%
drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
mutate(
Age = age_at_bs_ab_initiation,
"Female sex" = sex_0_male_1_female,
bsAb = bs_ab_name,
"Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
"ECOG" = case_when(
ecog_ps_at_bs_ab_initiation_0_5 %in% c("0","1") ~"0-1",
ecog_ps_at_bs_ab_initiation_0_5 ==2 ~"2",
ecog_ps_at_bs_ab_initiation_0_5 >=3 ~"3+",
.default = NULL),
EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
Hgb = as.double(baseline_hgb),
LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
"Prior BCMA" = prior_bcma_directed_therapy_yes_no
)
## Use the WeightIt function to generate propensity scores and weights
W.out <- weightit(bsAb ~ Age + `ECOG` + EMD + Hgb + LDH + `Prior BCMA`,
data=IPTW.data,
method = "glm", ## Methods: glm, gbm, cbps, npcbps, ebal, optweight, super, bart, energy
estimand = "ATE",
stabilize = TRUE)
bal.tab(W.out, stats = c("m", "v"), thresholds = c(m=0.25), binary = "std")
## Balance Measures
## Type Diff.Adj M.Threshold V.Ratio.Adj
## prop.score Distance 0.0050 Balanced, <0.25 0.9733
## Age Contin. 0.0090 Balanced, <0.25 0.9595
## ECOG_0-1 Binary -0.0121 Balanced, <0.25 .
## ECOG_2 Binary 0.0152 Balanced, <0.25 .
## ECOG_3+ Binary -0.0027 Balanced, <0.25 .
## EMD Binary -0.0034 Balanced, <0.25 .
## Hgb Contin. -0.0186 Balanced, <0.25 1.0974
## LDH Contin. -0.0276 Balanced, <0.25 1.0961
## Prior BCMA Binary -0.0119 Balanced, <0.25 .
##
## Balance tally for mean differences
## count
## Balanced, <0.25 9
## Not Balanced, >0.25 0
##
## Variable with the greatest mean difference
## Variable Diff.Adj M.Threshold
## LDH -0.0276 Balanced, <0.25
##
## Effective sample sizes
## elra tec
## Unadjusted 123. 248.
## Adjusted 113.38 242.39
summary(W.out)
## Summary of weights
##
## - Weight ranges:
##
## Min Max
## elra 0.4745 |---------------------------| 2.3134
## tec 0.7811 |-------------| 1.7418
##
## - Units with the 5 most extreme weights by group:
##
## 123 102 81 91 9
## elra 1.5458 1.5896 1.653 1.7607 2.3134
## 169 71 274 258 271
## tec 1.4257 1.4453 1.5215 1.6783 1.7418
##
## - Weight statistics:
##
## Coef of Var MAD Entropy # Zeros
## elra 0.292 0.222 0.040 0
## tec 0.152 0.115 0.011 0
##
## - Effective Sample Sizes:
##
## elra tec
## Unweighted 123. 248.
## Weighted 113.38 242.39
love_plot <- love.plot(W.out,
binary = "std",
line = TRUE,
abs = TRUE,
thresholds = c(m=0.2),
sample.names = c("Unweighted", "IPTW"),
title = "A)",
var.order = "unadjusted",
var.names = c(
prop.score = "Distance"
))+ theme(plot.title = element_text(hjust = 0))
love_plot
weight_df = data.frame(W.out$weights, W.out$treat)
IPTW.data$weights <- W.out$weights
IPTW.data$ps <- W.out$ps
write.csv(data1, "IPTW.dataset.csv")
Weights_plot <- ggplot(
IPTW.data, aes(x=ps, y = weights, color = bsAb)) +
geom_point(shape = 16) +
theme_classic() +
scale_y_continuous(breaks = seq(0,4 , by = 0.5)) +
labs(x = "Propensity score", y = "Weights") +
ggtitle("B)")
Weights_plot
theme_gtsummary_compact()
combined_db %>%
transmute(
"Site" = site,
bs_ab_name,
) %>%
tbl_summary(
by = bs_ab_name,
missing = "no",
sort = list(all_categorical() ~ "frequency"),
statistic = list(
all_continuous() ~ "{median} ({p25}, {p75})",
all_categorical() ~ "{n} ({p}%)"
)
) %>%
bold_labels()
| Characteristic | elra N = 1301 |
tec N = 2771 |
|---|---|---|
| Site | ||
| Â Â Â Â MCC | 53 (41%) | 72 (26%) |
| Â Â Â Â MDACC | 12 (9.2%) | 57 (21%) |
| Â Â Â Â CC | 20 (15%) | 47 (17%) |
| Â Â Â Â UTSW | 10 (7.7%) | 40 (14%) |
| Â Â Â Â FHCC | 7 (5.4%) | 24 (8.7%) |
| Â Â Â Â Huntsman | 9 (6.9%) | 22 (7.9%) |
| Â Â Â Â Kansas | 4 (3.1%) | 13 (4.7%) |
| Â Â Â Â HUMC | 10 (7.7%) | 1 (0.4%) |
| Â Â Â Â Roswell | 5 (3.8%) | 1 (0.4%) |
| 1 n (%) | ||
theme_gtsummary_compact()
data1 = IPTW.data %>%
transmute(
"Age (years)" = age_at_bs_ab_initiation,
"Female sex" = sex_0_male_1_female == 1,
"Race" = factor(case_when(
is.na(race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5) ~ NA,
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 0 ~ "White",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 1 ~ "Black",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 2 ~ "AAPI",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 3 ~ "American Indian",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 4 ~ "Other",
race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 5 ~ NA
), levels = c("White","Black","AAPI","American Indian","Other",NA)),
#ecog_ps_at_bs_ab_initiation_0_5,
"ECOG" = case_when(
ecog_ps_at_bs_ab_initiation_0_5 == 0 ~ "0",
ecog_ps_at_bs_ab_initiation_0_5 == 1 ~ "1",
ecog_ps_at_bs_ab_initiation_0_5 == 2 ~ "2",
ecog_ps_at_bs_ab_initiation_0_5 >= 3 ~ "3+"
),
"LDH (U/L)" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
"Hgb (g/dL)" = baseline_hgb,
"Baseline ferritin (ng/mL)" = as.double(baseline_ferritin_prior_to_bs_ab_initiation),
"Heavy chain" = factor(case_when(
myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgG" ~ "IgG",
myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgM" ~ "IgM",
myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgA" ~ "IgA",
myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgD" ~ "IgD",
.default = "None"
), levels = c("IgG", "IgA","IgM","IgD","None")),
"True EMD" = just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,
"≥1 HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes==1 | t_4_14_prior_to_bs_ab_0_no_1_yes == 1 | t_14_16_prior_to_bs_ab_0_no_1_yes == 1 | amp_1_q_only_0_no_1_yes == 1 | gain_amp1q21_prior_to_bs_ab_0_no_1_yes == 1,1,0),
"Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
"Prior ASCT" = prior_auto_sct_no_0_yes_1,
"Prior CAR-T" = prior_cart_0_no_1_yes,
"Triple refractory" = triple_class_refractory_i_mi_d_pi_anti_cd38 == 1,
"Penta refractory" = penta_refractory_2_p_is_2_i_mi_ds_anti_cd38 == 1,
"Prior BCMA-directed therapy" = !is.na(most_recent_bcma_agent_name),
"Trial eligible" = eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,
"Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
bs_ab_name,
weights
)
table_unweighted <- data1 %>% select(-weights) %>%
tbl_summary(
by = bs_ab_name,
missing = "no",
#sort = list(all_categorical() ~ "frequency"),
statistic = list(
all_continuous() ~ "{median} ({p25}, {p75})",
all_categorical() ~ "{n} ({p}%)",
"Age (years)" ~ "{median} ({min} to {max})"
)
) %>%
add_p() %>%
bold_labels()
table_IPTW <- survey::svydesign(
id = ~0,
weights = ~weights,
data = data1,
fpc = NULL
) %>%
tbl_svysummary(
by = bs_ab_name,
missing = "no",
statistic = list(
all_continuous() ~ "{median} ({p25}, {p75})",
all_categorical() ~ "{n} ({p}%)"),
include = c("Age (years)", "Female sex","Race","ECOG","LDH (U/L)","Hgb (g/dL)","Baseline ferritin (ng/mL)", "Heavy chain","True EMD","≥1 HRCA","Prior LOTs","Prior ASCT","Prior CAR-T","Triple refractory","Penta refractory","Trial eligible","Prior BCMA-directed therapy", "Received full dose")
) %>%
add_p() %>%
bold_labels()
tbl_merge(list(table_IPTW, table_unweighted), tab_spanner = c("**IPTW**", "**Unweighted**"))
| Characteristic |
IPTW
|
Unweighted
|
||||
|---|---|---|---|---|---|---|
| elra N = 1231 |
tec N = 2481 |
p-value2 | elra N = 1233 |
tec N = 2483 |
p-value4 | |
| Age (years) | 70 (64, 76) | 69 (62, 76) | 0.9 | 71 (39 to 95) | 68 (35 to 88) | 0.060 |
| Female sex | 66 (54%) | 114 (46%) | 0.2 | 68 (55%) | 115 (46%) | 0.11 |
| Race | 0.2 | 0.15 | ||||
| Â Â Â Â White | 94 (79%) | 187 (76%) | 95 (81%) | 187 (76%) | ||
| Â Â Â Â Black | 17 (14%) | 50 (20%) | 16 (14%) | 51 (21%) | ||
| Â Â Â Â AAPI | 3 (2.5%) | 6 (2.3%) | 3 (2.5%) | 5 (2.0%) | ||
| Â Â Â Â American Indian | 2 (1.8%) | 0 (0%) | 2 (1.7%) | 0 (0%) | ||
| Â Â Â Â Other | 2 (2.0%) | 4 (1.8%) | 2 (1.7%) | 4 (1.6%) | ||
| ECOG | >0.9 | 0.5 | ||||
| Â Â Â Â 0 | 23 (18%) | 49 (20%) | 21 (17%) | 50 (20%) | ||
| Â Â Â Â 1 | 62 (50%) | 119 (48%) | 59 (48%) | 121 (49%) | ||
| Â Â Â Â 2 | 28 (23%) | 59 (24%) | 29 (24%) | 60 (24%) | ||
| Â Â Â Â 3+ | 10 (8.2%) | 20 (8.2%) | 14 (11%) | 17 (6.9%) | ||
| LDH (U/L) | 223 (178, 280) | 208 (167, 290) | 0.4 | 222 (175, 277) | 208 (167, 291) | 0.5 |
| Hgb (g/dL) | 9.90 (8.50, 11.40) | 9.80 (8.40, 11.50) | 0.9 | 10.00 (8.70, 11.70) | 9.60 (8.25, 11.20) | 0.027 |
| Baseline ferritin (ng/mL) | 346 (140, 1,250) | 379 (125, 1,219) | >0.9 | 293 (137, 832) | 393 (128, 1,254) | 0.3 |
| Heavy chain | 0.001 | <0.001 | ||||
| Â Â Â Â IgG | 72 (58%) | 109 (44%) | 72 (59%) | 107 (43%) | ||
| Â Â Â Â IgA | 26 (21%) | 43 (17%) | 26 (21%) | 42 (17%) | ||
| Â Â Â Â IgM | 1 (0.7%) | 1 (0.4%) | 1 (0.8%) | 1 (0.4%) | ||
| Â Â Â Â IgD | 2 (1.5%) | 0 (0%) | 2 (1.6%) | 0 (0%) | ||
| Â Â Â Â None | 23 (19%) | 95 (38%) | 22 (18%) | 98 (40%) | ||
| True EMD | 20 (16%) | 39 (16%) | >0.9 | 26 (21%) | 34 (14%) | 0.067 |
| ≥1 HRCA | 81 (70%) | 151 (65%) | 0.3 | 81 (70%) | 152 (65%) | 0.4 |
| Prior LOTs | 6.00 (5.00, 8.00) | 6.00 (4.00, 8.00) | 0.4 | 6.00 (4.00, 8.00) | 6.00 (4.00, 8.00) | >0.9 |
| Prior ASCT | 66 (54%) | 141 (57%) | 0.6 | 61 (50%) | 144 (58%) | 0.12 |
| Prior CAR-T | 53 (43%) | 90 (36%) | 0.2 | 50 (41%) | 94 (38%) | 0.6 |
| Triple refractory | 112 (91%) | 224 (91%) | 0.9 | 112 (91%) | 226 (91%) | >0.9 |
| Penta refractory | 61 (49%) | 123 (49%) | >0.9 | 61 (50%) | 126 (51%) | 0.8 |
| Trial eligible | 24 (19%) | 55 (22%) | 0.5 | 26 (21%) | 53 (21%) | >0.9 |
| Prior BCMA-directed therapy | 67 (54%) | 133 (54%) | >0.9 | 60 (49%) | 138 (56%) | 0.2 |
| Received full dose | 114 (93%) | 241 (97%) | 0.059 | 114 (93%) | 241 (97%) | 0.045 |
| 1 Median (Q1, Q3); n (%) | ||||||
| 2 Design-based KruskalWallis test; Pearson’s X^2: Rao & Scott adjustment | ||||||
| 3 Median (Min to Max); n (%); Median (Q1, Q3) | ||||||
| 4 Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test | ||||||
data1 = IPTW.data %>%
transmute(
BsAb = as.factor(bs_ab_name),
CRS.any = ifelse( max_crs_grade_0_5 != 0, 1, 0),
ICANS.any = ifelse( max_icans_grade != 0,1,0),
CRS.two = ifelse( max_crs_grade_0_5 == 2 | max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
ICANS.two = ifelse( max_icans_grade ==2 | max_icans_grade == 3 | max_icans_grade == 4,1,0),
CRS.three = ifelse( max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
ICANS.three = ifelse( max_icans_grade == 3 | max_icans_grade == 4,1,0),
weights
)
weighted_table <- survey::svydesign(
id = ~0,
weights = ~weights,
data = data1,
fpc = NULL
)
combined_data <- rbind(
IPTW.data %>%
count(bs_ab_name, max_crs_grade_0_5, wt=weights) %>%
group_by(bs_ab_name) %>%
mutate(
bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
total = sum(n[max_crs_grade_0_5 != 0]),
pct = prop.table(n) * 100,
total_pct = sum(pct[max_crs_grade_0_5 != 0]),
Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
Category = "CRS"
),
IPTW.data %>%
count(bs_ab_name, max_icans_grade, wt=weights) %>%
group_by(bs_ab_name) %>%
mutate(
bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
total = sum(n[max_icans_grade != 0]),
pct = prop.table(n) * 100,
total_pct = sum(pct[max_icans_grade != 0]),
Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
Category = "ICANS"
)
) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS"))
)
# Any grade toxicity
CRS.pval <- as.numeric(svychisq(~CRS.any+BsAb, weighted_table)$p.value)
ICANS.pval <- as.numeric(svychisq(~ICANS.any+BsAb, weighted_table)$p.value)
f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS"))
)
# G2+ toxicity
CRS.two.pval <- as.numeric(svychisq(~CRS.two+BsAb, weighted_table)$p.value)
ICANS.two.pval <- as.numeric(svychisq(~ICANS.two+BsAb, weighted_table)$p.value)
f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS"))
)
# G3+ toxicity
CRS.three.pval <- as.numeric(svychisq(~CRS.three+BsAb, weighted_table)$p.value)
ICANS.three.pval <- as.numeric(svychisq(~ICANS.three+BsAb, weighted_table)$p.value)
f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS"))
)
toxplot_IPTW <- ggplot(combined_data) +
aes(x = bs_ab_name, y = pct, fill = Grade) +
geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")),
data = . %>% filter(Grade != 0 & Grade != 4),
position = position_stack(vjust = 0.5)) +
geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 6, vjust = 0),
data = . %>% filter(Grade == 1),
color = "#104a8e",
fontface = "bold") +
facet_wrap(~Category, scales = "free_y") +
ylab("Proportion of patients (% [n])") +
xlab(NULL) +
theme_classic()+
scale_fill_brewer(palette = "Pastel1") +
#guides(fill = guide_legend(reverse = TRUE)) +
scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
theme(
plot.title = element_text(size = 14),
axis.title = element_text(size = 14),
axis.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
strip.text = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1),
#legend.position = "bottom"
legend.position = "right"
) +
coord_cartesian(ylim = c(0, 115)) +
geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) +
geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) +
geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) +
ggtitle("A) IPTW")
toxplot_IPTW
combined_data <- rbind(
combined_db %>%
count(bs_ab_name, max_crs_grade_0_5) %>%
group_by(bs_ab_name) %>%
mutate(
total = sum(n[max_crs_grade_0_5 != 0]),
pct = prop.table(n) * 100,
total_pct = sum(pct[max_crs_grade_0_5 != 0]),
Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
Category = "CRS"
),
combined_db %>%
count(bs_ab_name, max_icans_grade) %>%
group_by(bs_ab_name) %>%
mutate(
total = sum(n[max_icans_grade != 0]),
pct = prop.table(n) * 100,
total_pct = sum(pct[max_icans_grade != 0]),
Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
Category = "ICANS"
)
) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS"))
)
## Any-grade toxicity
CRS.d.test <- data.frame(
results = ifelse(combined_db$max_crs_grade_0_5 !=0,1,0) ,
test = combined_db$bs_ab_name
)
CRS.pval <- chisq.test(CRS.d.test$results,CRS.d.test$test, correct = FALSE)$p.value
ICANS.d.test <- data.frame(
results = ifelse(combined_db$max_icans_grade !=0,1,0) ,
test = combined_db$bs_ab_name
)
ICANS.pval <- chisq.test(ICANS.d.test$results,ICANS.d.test$test, correct = FALSE)$p.value
f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS") )
)
## G2+ toxicity
CRS.two.d.test <- data.frame(
results = ifelse(combined_db$max_crs_grade_0_5 == 2 | combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
test = combined_db$bs_ab_name
)
CRS.two.pval <- chisq.test(CRS.two.d.test$results,CRS.two.d.test$test, correct = FALSE)$p.value
ICANS.two.d.test <- data.frame(
results = ifelse(combined_db$max_icans_grade == 2 | combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
test = combined_db$bs_ab_name
)
ICANS.two.pval <- chisq.test(ICANS.two.d.test$results,ICANS.two.d.test$test, correct = FALSE)$p.value
f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS") )
)
## G3+ toxicity
CRS.three.d.test <- data.frame(
results = ifelse( combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
test = combined_db$bs_ab_name
)
CRS.three.pval <- fisher.test(CRS.three.d.test$results,CRS.three.d.test$test)$p.value
ICANS.three.d.test <- data.frame(
results = ifelse(combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
test = combined_db$bs_ab_name
)
ICANS.three.pval <- chisq.test(ICANS.three.d.test$results,ICANS.three.d.test$test, correct = FALSE)$p.value
f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
mutate(
Category = factor(Category, levels = c("CRS","ICANS") )
)
toxplot_unweighted <- ggplot(combined_data) +
aes(x = bs_ab_name, y = pct, fill = Grade) +
geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")),
data = . %>% filter(Grade != 0 & Grade != 4),
position = position_stack(vjust = 0.5)) +
geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 6, vjust = 0),
data = . %>% filter(Grade == 1),
color = "#104a8e",
fontface = "bold") +
facet_wrap(~Category, scales = "free_y") +
ylab("Proportion of patients (% [n])") +
xlab(NULL) +
theme_classic()+
scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 100)) +
theme(
plot.title = element_text(size = 14),
axis.title = element_text(size = 14),
axis.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
strip.text = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none"
) +
coord_cartesian(ylim = c(0, 114)) +
scale_fill_brewer(palette = "Pastel1", direction = 1) +
geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) +
geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) +
geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) +
ggtitle("B) Unweighted")
toxplot_unweighted
FIGURE1 <- toxplot_IPTW / toxplot_unweighted
ggsave("svg/FIGURE1.svg", FIGURE1, device = "svg",
width = 7,
height = 11,
units = c("in"))
ggsave("tiff/FIGURE1.tiff", FIGURE1, device = "tiff",
width = 7,
height = 11,
dpi = 600,
units = c("in"),
compression = "lzw")
ggsave("jpeg/FIGURE1.jpeg", FIGURE1, device = "jpeg",
width = 7,
height = 11,
dpi = 600,
units = "in")
knitr::include_graphics("jpeg/FIGURE1.jpeg")
data1 = IPTW.data %>%
filter(!is.na(best_response)) %>%
transmute(
BsAb = as.factor(bs_ab_name),
best_response = factor(best_response, levels = c("sCR", "CR","VGPR","PR","MR","SD","PD")),
ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR",1,0),
VGPR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR",1,0),
CR = ifelse(best_response == "sCR" | best_response == "CR",1,0),
weights
)
weighted_table <- survey::svydesign(
id = ~0,
weights = ~weights,
data = data1,
fpc = NULL
)
data1 <- IPTW.data %>%
filter(!is.na(best_response)) %>%
count(bs_ab_name, best_response, wt=weights) %>%
group_by(bs_ab_name) %>%
mutate(
total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
pct = prop.table(n) * 100,
total_pct = sum(pct[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
best_response = factor(best_response, levels = c("sCR", "CR","VGPR","PR","MR", "SD","PD")),
bs_ab_name = as.factor(bs_ab_name),
Category = "Response rate"
)
Response.pval <- as.numeric(svychisq(~ORR+BsAb, weighted_table)$p.value)
CR.pval <- as.numeric(svychisq(~CR+BsAb, weighted_table)$p.value)
VGPR.pval <- as.numeric(svychisq(~VGPR+BsAb, weighted_table)$p.value)
f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
mutate(
Category = factor(Category, levels = c("Response rate"))
)
f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
mutate(
Category = factor(Category, levels = c("Response rate"))
)
f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
mutate(
Category = factor(Category, levels = c("Response rate"))
)
responseplot_IPTW <- ggplot(data1) +
aes(x = bs_ab_name, y = pct, fill = best_response) +
geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR")) +
geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")),
data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR"| best_response == "PR"),
position = position_stack(vjust = 0.5)) +
geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 7, vjust = 0),
data = . %>% filter(best_response == "PR"),
color = "#104a8e",
fontface = "bold") +
facet_wrap(~Category, scales = "free_y") +
ylab("Proportion of patients (% [n])") +
xlab(NULL) +
theme_classic()+
scale_fill_brewer(
palette = "Pastel1",
direction = -1,
breaks = rev(levels(data1$best_response))
) +
guides(fill = guide_legend(reverse = TRUE)) +
scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
theme(
plot.title = element_text(size = 14),
axis.title = element_text(size = 14),
axis.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
strip.text = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1),
#legend.position = "bottom"
legend.position = "right"
) +
coord_cartesian(ylim = c(0, 114)) +
geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) +
geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR ) +
geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR ) +
labs(fill = "Best response")+
ggtitle("A) IPTW")
responseplot_IPTW
data1 <- combined_db %>%
filter(!is.na(best_response)) %>%
transmute(
best_response = factor(
best_response,
levels = c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" |best_response == "PR",1,0),
bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
Category = "Response rate"
)
combined_data <- combined_db %>%
filter(!is.na(best_response)) %>%
count(bs_ab_name, best_response) %>%
group_by(bs_ab_name) %>%
transmute(
bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
best_response = factor(
best_response,
levels = c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
pct = prop.table(n) * 100,
n,
total_pct = sum(pct[best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR"]),
Category = "Response rate"
)
## ORR
Response.d.test <- data.frame(
results = ifelse(data1$ORR,1,0),
test = data1$bs_ab_name
)
Response.pval <- chisq.test(Response.d.test$results,Response.d.test$test, correct = FALSE)$p.value
## VGPR
VGPR.d.test <- data.frame(
results = ifelse(data1$best_response == "CR" | data1$best_response == "VGPR",1,0),
test = data1$bs_ab_name
)
VGPR.pval <- chisq.test(VGPR.d.test$results,VGPR.d.test$test, correct = FALSE)$p.value
## CR
CR.d.test <- data.frame(
results = ifelse(data1$best_response == "CR",1,0),
test = data1$bs_ab_name
)
CR.pval <- chisq.test(CR.d.test$results,CR.d.test$test, correct = FALSE)$p.value
f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
mutate(
Category = factor(Category, levels = c("Response rate"))
)
f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
mutate(
Category = factor(Category, levels = c("Response rate"))
)
f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
mutate(
Category = factor(Category, levels = c("Response rate"))
)
responseplot_unweighted <- ggplot(combined_data) +
aes(x = bs_ab_name, y = pct, fill = best_response) +
geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR")) +
geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")),
data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"),
position = position_stack(vjust = 0.5)) +
geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 7, vjust = 0),
data = . %>% filter(best_response == "PR"),
color = "#104a8e",
fontface = "bold") +
facet_wrap(~Category, scales = "free_y") +
ylab("Proportion of patients (% [n])") +
xlab(NULL) +
theme_classic()+
scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
theme(
plot.title = element_text(size = 14),
axis.title = element_text(size = 14),
axis.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
strip.text = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1),
#legend.position = "bottom"
legend.position = "none"
) +
coord_cartesian(ylim = c(0, 114)) +
scale_fill_brewer(palette = "Pastel1", direction = -1) +
geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels) +
geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR ) +
geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR ) +
labs(fill = "Best response")+
ggtitle("B) Unweighted")
responseplot_unweighted
FIGURE2 <- responseplot_IPTW / responseplot_unweighted
ggsave("svg/FIGURE2.svg", FIGURE2, device = "svg",
width = 7,
height = 11,
units = c("in"))
ggsave("tiff/FIGURE2.tiff", FIGURE2, device = "tiff",
width = 7,
height = 11,
dpi = 600,
units = c("in"),
compression = "lzw")
ggsave("jpeg/FIGURE2.jpeg", FIGURE2, device = "jpeg",
width = 7,
height = 11,
dpi = 600,
units = "in")
knitr::include_graphics("jpeg/FIGURE2.jpeg")
data1 <- IPTW.data %>%
transmute(
bs_ab_name,
Death = ifelse(death_yes_no,1,0),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Reverse_death = ifelse(death_yes_no == 1, 0,1),
weights
)
reverse_km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1)
reverse_km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~
## 1, data = data1, weights = weights)
##
## 15 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## [1,] 356 356 215 13.4 12.5 14.9
reverse_km_OS_elra <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "elra"))
reverse_km_OS_elra
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~
## 1, data = data1 %>% filter(bs_ab_name == "elra"), weights = weights)
##
## 6 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## [1,] 117 117 74.5 7.47 6.03 9
reverse_km_OS_tec <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "tec"))
reverse_km_OS_tec
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~
## 1, data = data1 %>% filter(bs_ab_name == "tec"), weights = weights)
##
## 9 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## [1,] 239 239 140 17.1 14.9 20.1
data1 <- IPTW.data %>%
transmute(
site,
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
bs_ab_name,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_bs_ab_step_up_d1,
date_of_pd,
date_of_last_contact,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
weights
)
data_dor <- IPTW.data %>%
filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
transmute(
weights,
bs_ab_name,
day_30_response,
day_90_response,
best_response,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_pd,
date_of_last_contact,
date_of_bs_ab_step_up_d1,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
day_30_response,
day_90_response,
best_response,
time_to_best_response = case_when(
best_response == day_30_response ~ 30,
best_response == day_90_response ~ 90,
.default = NA
),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
) %>%
filter(!is.na(time_to_best_response)) %>%
mutate(
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1)- time_to_best_response,
) %>%
filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)
km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name,
## data = data1, weights = weights)
##
## 15 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 117 117 42.8 NA 8.6 NA
## bs_ab_name=tec 239 239 98.6 21.5 16.4 NA
med_time <- surv_median(km_OS)
cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))
OS <- ggsurvplot(km_OS,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Overall survival",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "A)",
size = 0.5
)
OS$plot <- OS$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
OS
km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~
## bs_ab_name, data = data1, weights = weights)
##
## 12 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 118 118 68.1 3.77 2.57 7.17
## bs_ab_name=tec 241 241 151.1 8.50 7.10 11.30
med_time <- surv_median(km_PFS)
cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))
PFS <- ggsurvplot(km_PFS,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Progression-free survival",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "B)",
size = 0.5
)
PFS$plot <- PFS$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
PFS
km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~
## bs_ab_name, data = data_dor, weights = weights)
##
## 3 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 62 60.2 20.6 11.7 6.13 NA
## bs_ab_name=tec 118 120.7 56.1 13.7 9.27 22.3
med_time <- surv_median(km_DOR)
cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))
DOR <- ggsurvplot(km_DOR,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
#risk.table.fontsize = 4,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Duration of response",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "C)",
size = 0.5
)
DOR$plot <- DOR$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
DOR
FIGURE3 <- wrap_plots(
(OS$plot / OS$table) + plot_layout(heights = c(4, 0.9)),
(PFS$plot / PFS$table) + plot_layout(heights = c(4, 0.9)),
(DOR$plot / DOR$table) + plot_layout(heights = c(4, 0.9)),
ncol = 1 # Force vertical stacking
)
ggsave("svg/FIGURE3.svg", FIGURE3, device = "svg",
width = 6.5,
height = 8.5,
units = c("in"))
ggsave("jpeg/FIGURE3.jpeg", FIGURE3, device = "jpeg",
width = 6.5,
height = 8.5,
dpi = 600,
units = c("in"))
ggsave("tiff/FIGURE3.tiff", FIGURE3, device = "tiff",
width = 6.5,
height = 8.5,
dpi = 600,
units = c("in"),
compression = "lzw")
knitr::include_graphics("jpeg/FIGURE3.jpeg")
data1 <- IPTW.data %>%
transmute(
site,
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
elra_ind = ifelse(bs_ab_name == "elra", 1, 0), # 1 = elra, 0 = tec
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_bs_ab_step_up_d1,
date_of_pd,
date_of_last_contact,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
weights
)
## ---- Cox models ----
cox_OS <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ elra_ind, weights = weights, data = data1)
cox_PFS <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ elra_ind, weights = weights, data = data1)
## ---- PH diagnostics ----
zph_OS <- cox.zph(cox_OS)
zph_PFS <- cox.zph(cox_PFS)
# Global PH p-values (same as elra_ind right now)
p_OS <- signif(zph_OS$table["GLOBAL", "p"], 3)
p_PFS <- signif(zph_PFS$table["GLOBAL", "p"], 3)
par(mfrow = c(1, 2))
plot(
zph_OS,
main = paste0(
"Schoenfeld residuals – OS\n",
"Global PH test: p = ", p_OS
)
)
plot(
zph_PFS,
main = paste0(
"Schoenfeld residuals – PFS\n",
"Global PH test: p = ", p_PFS
)
)
par(mfrow = c(1, 1))
data1 <- IPTW.data %>%
filter(EMD == 1) %>%
transmute(
site,
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
bs_ab_name,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_bs_ab_step_up_d1,
date_of_pd,
date_of_last_contact,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
weights
)
data_dor <- IPTW.data %>%
filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
filter(EMD == 1) %>%
transmute(
weights,
bs_ab_name,
day_30_response,
day_90_response,
best_response,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_pd,
date_of_last_contact,
date_of_bs_ab_step_up_d1,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
day_30_response,
day_90_response,
best_response,
time_to_best_response = case_when(
best_response == day_30_response ~ 30,
best_response == day_90_response ~ 90,
.default = NA
),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
) %>%
filter(!is.na(time_to_best_response)) %>%
mutate(
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1)- time_to_best_response,
) %>%
filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)
km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name,
## data = data1, weights = weights)
##
## 4 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 24 18.2 4.86 NA 5.23 NA
## bs_ab_name=tec 32 37.4 21.05 14.1 10.43 NA
med_time <- surv_median(km_OS)
cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))
OS <- ggsurvplot(km_OS,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Overall survival",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "A)",
size = 0.5
)
OS$plot <- OS$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
OS
km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~
## bs_ab_name, data = data1, weights = weights)
##
## 1 observation deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 25 18.8 8.27 6.27 2.57 NA
## bs_ab_name=tec 34 39.5 27.45 7.93 4.67 12.5
med_time <- surv_median(km_PFS)
cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))
PFS <- ggsurvplot(km_PFS,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Progression-free survival",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "B)",
size = 0.5
)
PFS$plot <- PFS$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
PFS
km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~
## bs_ab_name, data = data_dor, weights = weights)
##
## 1 observation deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 11 8.76 2.39 NA 5.27 NA
## bs_ab_name=tec 15 18.05 8.69 6.93 5.63 NA
med_time <- surv_median(km_DOR)
cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))
DOR <- ggsurvplot(km_DOR,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
#risk.table.fontsize = 4,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Duration of response",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "C)",
size = 0.5
)
DOR$plot <- DOR$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
DOR
data1 <- IPTW.data %>%
filter(ecog_ps_at_bs_ab_initiation_0_5 >= 2) %>%
transmute(
site,
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
bs_ab_name,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_bs_ab_step_up_d1,
date_of_pd,
date_of_last_contact,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
weights
)
data_dor <- IPTW.data %>%
filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
filter(ecog_ps_at_bs_ab_initiation_0_5 >= 2) %>%
transmute(
weights,
bs_ab_name,
day_30_response,
day_90_response,
best_response,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_pd,
date_of_last_contact,
date_of_bs_ab_step_up_d1,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
day_30_response,
day_90_response,
best_response,
time_to_best_response = case_when(
best_response == day_30_response ~ 30,
best_response == day_90_response ~ 90,
.default = NA
),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
) %>%
filter(!is.na(time_to_best_response)) %>%
mutate(
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1)- time_to_best_response,
) %>%
filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)
km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name,
## data = data1, weights = weights)
##
## 4 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 42 37.9 18.8 4.87 2.27 NA
## bs_ab_name=tec 74 75.8 45.8 9.57 5.27 18.8
med_time <- surv_median(km_OS)
cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))
OS <- ggsurvplot(km_OS,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Overall survival",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "A)",
size = 0.5
)
OS$plot <- OS$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
OS
km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~
## bs_ab_name, data = data1, weights = weights)
##
## 2 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 43 38.6 27.9 1.63 1.07 10.27
## bs_ab_name=tec 75 76.7 56.7 4.43 2.10 8.73
med_time <- surv_median(km_PFS)
cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))
PFS <- ggsurvplot(km_PFS,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Progression-free survival",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "B)",
size = 0.5
)
PFS$plot <- PFS$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
PFS
km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~
## bs_ab_name, data = data_dor, weights = weights)
##
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 15 12.9 6.02 11.7 1.0 NA
## bs_ab_name=tec 30 31.3 16.87 13.3 8.3 NA
med_time <- surv_median(km_DOR)
cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))
DOR <- ggsurvplot(km_DOR,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
#risk.table.fontsize = 4,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Duration of response",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "C)",
size = 0.5
)
DOR$plot <- DOR$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
DOR
data1 <- IPTW.data %>%
filter(baseline_hgb <10) %>%
transmute(
site,
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
bs_ab_name,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_bs_ab_step_up_d1,
date_of_pd,
date_of_last_contact,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
weights
)
data_dor <- IPTW.data %>%
filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
filter(baseline_hgb <10) %>%
transmute(
weights,
bs_ab_name,
day_30_response,
day_90_response,
best_response,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_pd,
date_of_last_contact,
date_of_bs_ab_step_up_d1,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
day_30_response,
day_90_response,
best_response,
time_to_best_response = case_when(
best_response == day_30_response ~ 30,
best_response == day_90_response ~ 90,
.default = NA
),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
) %>%
filter(!is.na(time_to_best_response)) %>%
mutate(
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1)- time_to_best_response,
) %>%
filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)
km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name,
## data = data1, weights = weights)
##
## 8 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 55 63.3 31.7 4.87 3.17 NA
## bs_ab_name=tec 134 125.3 71.5 12.83 9.57 18.6
med_time <- surv_median(km_OS)
cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))
OS <- ggsurvplot(km_OS,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Overall survival",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "A)",
size = 0.5
)
OS$plot <- OS$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
OS
km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~
## bs_ab_name, data = data1, weights = weights)
##
## 4 observations deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 56 63.9 43.7 1.63 1.10 6.27
## bs_ab_name=tec 137 128.3 97.8 3.70 2.17 7.40
med_time <- surv_median(km_PFS)
cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))
PFS <- ggsurvplot(km_PFS,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Progression-free survival",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "B)",
size = 0.5
)
PFS$plot <- PFS$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
PFS
km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~
## bs_ab_name, data = data_dor, weights = weights)
##
## 1 observation deleted due to missingness
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 26 28.6 13.6 6.17 4.2 NA
## bs_ab_name=tec 51 47.5 27.4 9.90 7.0 NA
med_time <- surv_median(km_DOR)
cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))
DOR <- ggsurvplot(km_DOR,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
#risk.table.fontsize = 4,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Duration of response",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "C)",
size = 0.5
)
DOR$plot <- DOR$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
DOR
data1 <- IPTW.data %>%
filter(age_at_bs_ab_initiation >=80) %>%
transmute(
site,
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
bs_ab_name,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_bs_ab_step_up_d1,
date_of_pd,
date_of_last_contact,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
weights
)
data_dor <- IPTW.data %>%
filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
filter(age_at_bs_ab_initiation >=80) %>%
transmute(
weights,
bs_ab_name,
day_30_response,
day_90_response,
best_response,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_pd,
date_of_last_contact,
date_of_bs_ab_step_up_d1,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
day_30_response,
day_90_response,
best_response,
time_to_best_response = case_when(
best_response == day_30_response ~ 30,
best_response == day_90_response ~ 90,
.default = NA
),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
) %>%
filter(!is.na(time_to_best_response)) %>%
mutate(
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1)- time_to_best_response,
) %>%
filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)
km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name,
## data = data1, weights = weights)
##
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 18 15.0 11.1 5.2 3.47 NA
## bs_ab_name=tec 30 32.9 10.6 17.7 9.53 NA
med_time <- surv_median(km_OS)
cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))
OS <- ggsurvplot(km_OS,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Overall survival",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "A)",
size = 0.5
)
OS$plot <- OS$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
OS
km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~
## bs_ab_name, data = data1, weights = weights)
##
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 18 15.0 11.1 3.73 1.63 NA
## bs_ab_name=tec 30 32.9 14.8 8.63 6.07 NA
med_time <- surv_median(km_PFS)
cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))
PFS <- ggsurvplot(km_PFS,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Progression-free survival",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "B)",
size = 0.5
)
PFS$plot <- PFS$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
PFS
km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~
## bs_ab_name, data = data_dor, weights = weights)
##
## records n events median 0.95LCL 0.95UCL
## bs_ab_name=elra 9 8.13 5.47 4.20 1.00 NA
## bs_ab_name=tec 19 20.84 7.68 6.63 5.07 NA
med_time <- surv_median(km_DOR)
cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)
HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)
## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))
DOR <- ggsurvplot(km_DOR,
pval=FALSE,
conf.int = FALSE,
risk.table = TRUE, # Add risk table
fontsize = 4,
risk.table.col = "strata", # Change risk table color by groups
tables.height = 0.3,
#risk.table.fontsize = 4,
tables.theme = theme_cleantable() +
theme(
axis.text.y = element_text(size = 9),
plot.title = element_text(size = 10)
),
xlab="Time (months)", ylab = "Duration of response",
xlim = c(0,18), break.x.by = c(3),
ylim = c(0,1), break.y.by = c(0.25),
legend = "none",
surv.scale="percent",
font.main = c(10, "plain", "black"),
font.x = c(10, "plain", "black"),
font.y = c(10, "plain", "black"),
font.caption = c(10, "plain", "black"),
font.tickslab = c(10, "plain", "black"),
legend.labs = c("Elra", "Tec"),
title = "C)",
size = 0.5
)
DOR$plot <- DOR$plot +
geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
annotate("text", x = -0.7, y = 0,
label = paste0("Elra vs tec\n",
" Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
" HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
color = "black", hjust = 0, vjust = 0, size = 3)
DOR
FIGURE4 <- wrap_plots(
(OS$plot / OS$table) + plot_layout(heights = c(4, 0.9)),
(PFS$plot / PFS$table) + plot_layout(heights = c(4, 0.9)),
(DOR$plot / DOR$table) + plot_layout(heights = c(4, 0.9)),
ncol = 1 # Force vertical stacking
)
ggsave("svg/FIGURE4.svg", FIGURE4, device = "svg",
width = 6.5,
height = 8.5,
units = c("in"))
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).
ggsave("jpeg/FIGURE4.jpeg", FIGURE4, device = "jpeg",
width = 6.5,
height = 8.5,
dpi = 600,
units = c("in"))
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).
ggsave("tiff/FIGURE4.tiff", FIGURE4, device = "tiff",
width = 6.5,
height = 8.5,
dpi = 600,
units = c("in"),
compression = "lzw")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).
knitr::include_graphics("jpeg/FIGURE4.jpeg")
options(gtsummary.add_p_footnote = FALSE, gtsummary.add_estimate_to_reference_rows = FALSE)
theme_gtsummary_compact()
preds0 <- c("Response", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds1 <- c("CR/sCR", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds2 <- c("BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+","True EMD", "Prior LOTs", "Prior BCMA")
preds3 <- c("ICANS 1+", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds4 <- c("CRS 1+", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds5 <- c("Steroid use", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds6 <- c("VGPR", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
data0 <- IPTW.data %>%
transmute(
weights,
BsAb = factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra"), levels = c("Tec","Elra")),
HT_Risk,
best_response,
"Infection" = ifelse(infection == 1, 1,0),
"CrCl <40" = ifelse(cr_cl_prior_to_step_up_dose_number_1 <40, 1,0),
"CR/sCR" = ifelse(best_response %in% c("CR","sCR"),1,0),
"VGPR" = ifelse(best_response %in% c("CR","sCR", "VGPR"),1,0),
"Response" = ifelse(best_response %in% c("CR","sCR", "VGPR", "PR"),1,0),
"True EMD" = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2, 1,0),
"Steroid use" = ifelse(steroid_use_as_tx==1,1,0),
"Age" = age_at_bs_ab_initiation/10,
"ECOG" = ecog_ps_at_bs_ab_initiation_0_5,
"ECOG 2+" = case_when(
is.na(ecog_ps_at_bs_ab_initiation_0_5) ~ NA,
ecog_ps_at_bs_ab_initiation_0_5 %in% c(2,3,4) ~ 1,
TRUE ~ 0),
"ICANS 1+" = case_when(
is.na(max_icans_grade) ~ NA,
max_icans_grade %in% c(1,2,3,4) ~1,
TRUE ~ 0
),
"CRS 1+" = case_when(
is.na(max_crs_grade_0_5) ~ NA,
max_crs_grade_0_5 %in% c(1,2,3,4) ~1,
TRUE ~ 0
),
"Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
"HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes | t_4_14_prior_to_bs_ab_0_no_1_yes | t_14_16_prior_to_bs_ab_0_no_1_yes | gain_amp1q21_prior_to_bs_ab_0_no_1_yes, 1, 0),
"Prior BCMA" = prior_bcma_directed_therapy_yes_no,
prior_bcma_directed_therapy_yes_no,
Days.since.last.BCMA = case_when(
is.na(prior_bcma_directed_therapy_yes_no) ~ NA,
date_of_last_exposure_to_prior_bcma_agent > date_of_bs_ab_step_up_d1 ~ NA,
prior_bcma_directed_therapy_yes_no == 1 ~ as.numeric(as.Date(date_of_bs_ab_step_up_d1) - as.Date(date_of_last_exposure_to_prior_bcma_agent))
),
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_pd,
date_of_last_contact,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1),
"ALC BL" = baseline_alc,
"CrCl BL" = cr_cl_prior_to_step_up_dose_number_1,
"Trial ineligible" = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,0,1),
"Baseline platelets" = baseline_pl_ts,
"Baseline Hgb" = baseline_hgb,
"Baseline LDH" = as.double(baseline_ldh_prior_to_bs_ab_initiation)/100,
"Baseline ferritin" = as.double(baseline_ldh_prior_to_bs_ab_initiation)/1000,
time_to_best_response = case_when(
best_response == day_30_response ~ 30,
best_response == day_90_response ~ 90,
.default = NA
),
date_of_bs_ab_step_up_d1,
bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes,
days_to_infection_or_last_contact = as.numeric(
pmin(as.Date(date_of_infection_74), as.Date(date_of_last_contact), na.rm = TRUE) - as.Date(date_of_bs_ab_step_up_d1)
),
infection_or_death = ifelse(!is.na(date_of_infection_74) | death_yes_no == 1,1,0),
IVIg = ifelse(ivig_use_yes_1_no_0 == 1, 1, 0),
ALPS = ifelse(baseline_hgb >10 & baseline_ldh_prior_to_bs_ab_initiation <300,1,0),
baseline_pl_ts,
baseline_anc,
baseline_hgb,
baseline_crp_prior_to_bs_ab_initiation,
baseline_ferritin_prior_to_bs_ab_initiation,
CAR.HEMATOTOX = baseline_pl_ts * -0.05 + baseline_anc * -0.03 + baseline_hgb * -0.20 + as.numeric(baseline_crp_prior_to_bs_ab_initiation) * 0.15 + as.numeric(baseline_ferritin_prior_to_bs_ab_initiation) * 0.001
) %>%
mutate(
"BCMA exposure status" = factor(case_when(
#prior_bcma_directed_therapy_yes_no == 0 ~ " Never",
Days.since.last.BCMA < 365 ~ " <1 yr",
Days.since.last.BCMA >= 365 ~ " ≥1 yr"
), c(" ≥1 yr"," <1 yr")),
`CAR-HEMATOTOX risk` = factor(
case_when(
CAR.HEMATOTOX < 0 ~ "Low",
CAR.HEMATOTOX <= 5 ~ "Intermediate",
CAR.HEMATOTOX > 5 ~ "High"
),
levels = c("Low", "Intermediate", "High")
)
)
data1 <- data0 %>% filter(bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes==0) ## Exclude patients who received bsAb as bridging
data2 <- data1 %>%
filter(best_response %in% c("PR","VGPR","sCR", "CR")) %>%
filter(!is.na(time_to_best_response)) %>%
mutate(
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1)- time_to_best_response,
) %>%
filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)
## Response
uv_tab <- tbl_uvregression(
data0 %>%
dplyr::select(all_of(preds0)),
method = glm,
y = `Response`,
method.args = list(family = binomial, weights = data0$weights),
exponentiate = TRUE,
add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")
mv_tab<-glm( `Response` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
dplyr::select(all_of(preds0), weights = weights) , family = binomial ) %>%
tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")
table_response <- tbl_merge(
list(uv_tab, mv_tab),
tab_spanner = c("**Univariate**", "**Multivariable**")
)
## VGPR or better
uv_tab <- tbl_uvregression(
data0 %>%
dplyr::select(all_of(preds6)),
method = glm,
y = `VGPR`,
method.args = list(family = binomial, weights = data0$weights),
exponentiate = TRUE,
add_estimate_to_reference_rows = FALSE,
) %>% modify_header(estimate = "**Estimate**")
mv_tab<-glm( `VGPR` ~ BsAb +`True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
dplyr::select(all_of(preds6), weights = weights) , family = binomial ) %>%
tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")
table_VGPR <- tbl_merge(
list(uv_tab, mv_tab),
tab_spanner = c("**Univariate**", "**Multivariable**")
)
## sCR/CR
uv_tab <- tbl_uvregression(
data0 %>%
dplyr::select(all_of(preds1)),
method = glm,
y = `CR/sCR`,
method.args = list(family = binomial, weights = data0$weights),
exponentiate = TRUE,
add_estimate_to_reference_rows = FALSE,
) %>% modify_header(estimate = "**Estimate**")
mv_tab<-glm( `CR/sCR` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
dplyr::select(all_of(preds1), weights = weights) , family = binomial ) %>%
tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")
table_sCR.CR <- tbl_merge(
list(uv_tab, mv_tab),
tab_spanner = c("**Univariate**", "**Multivariable**")
)
## OS
uv_tab <- tbl_uvregression(
data1[c(preds2)],
method = coxph,
y = Surv(data1$Days.to.death.or.DLC, data1$Death),
exponentiate = TRUE,
method.args = list(weights = data1$weights)
) %>% modify_header(estimate = "**Estimate**")
mv_tab<-coxph( Surv(Days.to.death.or.DLC, Death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data1, weights = weights) %>%
tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")
table_os <- tbl_merge(
list(uv_tab, mv_tab),
tab_spanner = c("**Univariate**", "**Multivariable**")
)
## PFS
uv_tab <- tbl_uvregression(
data1[c(preds2)],
method = coxph,
y = Surv(data1$Days.to.death.or.relapse.or.DLC, data1$Relapse.or.death),
exponentiate = TRUE,
method.args = list(weights = data1$weights)
) %>% modify_header(estimate = "**Estimate**")
mv_tab<-coxph( Surv(Days.to.death.or.relapse.or.DLC, Relapse.or.death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data1) %>%
tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")
table_pfs <- tbl_merge(
list(uv_tab, mv_tab),
tab_spanner = c("**Univariate**", "**Multivariable**")
)
## DOR
uv_tab <- tbl_uvregression(
data2[c(preds2)],
method = coxph,
y = Surv(data2$Days.to.death.or.relapse.or.DLC, data2$Relapse.or.death),
exponentiate = TRUE,
method.args = list(weights = data2$weights)
) %>% modify_header(estimate = "**Estimate**")
mv_tab<-coxph( Surv( Days.to.death.or.relapse.or.DLC, Relapse.or.death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data2) %>%
tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")
table_dor <- tbl_merge(
list(uv_tab, mv_tab),
tab_spanner = c("**Univariate**", "**Multivariable**")
)
## IFS
uv_tab <- tbl_uvregression(
data1[c(preds2)],
method = coxph,
y = Surv(data1$days_to_infection_or_last_contact, data1$infection_or_death),
exponentiate = TRUE,
method.args = list(weights = data1$weights)
) %>% modify_header(estimate = "**Estimate**")
mv_tab<-coxph( Surv( days_to_infection_or_last_contact, infection_or_death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data1) %>%
tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")
table_ifs <- tbl_merge(
list(uv_tab, mv_tab),
tab_spanner = c("**Univariate**", "**Multivariable**")
)
## ICANS 1+
uv_tab <- tbl_uvregression(
data0 %>%
dplyr::select(all_of(preds3)),
method = glm,
y = `ICANS 1+`,
method.args = list(family = binomial, weights = data0$weights),
exponentiate = TRUE,
add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")
mv_tab<-glm( `ICANS 1+` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
dplyr::select(all_of(preds3), weights = weights) , family = binomial ) %>%
tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")
table_ICANS <- tbl_merge(
list(uv_tab, mv_tab),
tab_spanner = c("**Univariate**", "**Multivariable**")
)
## CRS 1+
uv_tab <- tbl_uvregression(
data0 %>%
dplyr::select(all_of(preds4)),
method = glm,
y = `CRS 1+`,
method.args = list(family = binomial, weights = data0$weights),
exponentiate = TRUE,
add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")
mv_tab<-glm( `CRS 1+` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
dplyr::select(all_of(preds4), weights = weights) , family = binomial ) %>%
tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")
table_CRS <- tbl_merge(
list(uv_tab, mv_tab),
tab_spanner = c("**Univariate**", "**Multivariable**")
)
## Steroid use
uv_tab <- tbl_uvregression(
data0 %>%
dplyr::select(all_of(preds5)),
method = glm,
y = `Steroid use`,
method.args = list(family = binomial, weights = data0$weights),
exponentiate = TRUE,
add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")
mv_tab<-glm( `Steroid use` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
dplyr::select(all_of(preds5), weights = weights) , family = binomial ) %>%
tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")
table_steroid <- tbl_merge(
list(uv_tab, mv_tab),
tab_spanner = c("**Univariate**", "**Multivariable**")
)
# Combine all tables into one
final_table <- tbl_stack(
list(
table_response %>% modify_header() %>% modify_footnote(everything() ~ NA, abbreviation = TRUE) %>%
modify_footnote(estimate_1 ~ "Estimate = Odds ratio or hazard ratio", abbreviation = TRUE),
table_VGPR,
table_sCR.CR,
table_os,
table_pfs,
table_dor,
table_ifs,
table_ICANS,
table_CRS,
table_steroid
),
group_header = c(
"Response",
"≥VGPR",
"CR/sCR",
"Overall Survival",
"Progression-free Survival",
"Duration of Response",
"Infection-free Survival",
"ICANS 1+",
"CRS 1+",
"Steroid use"
)
)
final_table
| Characteristic |
Univariate
|
Multivariable
|
|||||
|---|---|---|---|---|---|---|---|
| N | Estimate1 | 95% CI | p-value | Estimate | 95% CI | p-value | |
| Response | |||||||
| BsAb | 370 | ||||||
|     Tec | — | — | — | — | |||
| Â Â Â Â Elra | 0.86 | 0.55, 1.35 | 0.5 | 0.86 | 0.52, 1.42 | 0.6 | |
| Baseline Hgb | 370 | 1.46 | 1.29, 1.65 | <0.001 | 1.34 | 1.18, 1.53 | <0.001 |
| Baseline LDH | 370 | 0.86 | 0.77, 0.94 | 0.002 | 0.86 | 0.77, 0.95 | 0.004 |
| ECOG 2+ | 370 | 0.39 | 0.25, 0.61 | <0.001 | 0.42 | 0.26, 0.69 | <0.001 |
| True EMD | 370 | 0.91 | 0.52, 1.63 | 0.8 | 0.94 | 0.50, 1.78 | 0.9 |
| Prior LOTs | 370 | 0.87 | 0.79, 0.94 | 0.001 | 0.91 | 0.82, 1.01 | 0.087 |
| Prior BCMA | 370 | 0.56 | 0.36, 0.85 | 0.008 | 0.50 | 0.29, 0.87 | 0.014 |
| ≥VGPR | |||||||
| BsAb | 370 | ||||||
|     Tec | — | — | — | — | |||
| Â Â Â Â Elra | 0.62 | 0.40, 0.95 | 0.031 | 0.57 | 0.35, 0.93 | 0.026 | |
| Baseline Hgb | 370 | 1.38 | 1.23, 1.55 | <0.001 | 1.32 | 1.17, 1.49 | <0.001 |
| Baseline LDH | 370 | 0.84 | 0.73, 0.93 | 0.003 | 0.84 | 0.73, 0.93 | 0.004 |
| ECOG 2+ | 370 | 0.52 | 0.33, 0.81 | 0.004 | 0.57 | 0.34, 0.93 | 0.026 |
| True EMD | 370 | 0.63 | 0.35, 1.10 | 0.11 | 0.56 | 0.30, 1.04 | 0.071 |
| Prior LOTs | 370 | 0.88 | 0.80, 0.95 | 0.003 | 0.95 | 0.86, 1.06 | 0.4 |
| Prior BCMA | 370 | 0.47 | 0.31, 0.72 | <0.001 | 0.41 | 0.24, 0.70 | 0.001 |
| CR/sCR | |||||||
| BsAb | 370 | ||||||
|     Tec | — | — | — | — | |||
| Â Â Â Â Elra | 0.70 | 0.44, 1.10 | 0.12 | 0.68 | 0.40, 1.12 | 0.13 | |
| Baseline Hgb | 370 | 1.44 | 1.29, 1.63 | <0.001 | 1.39 | 1.23, 1.58 | <0.001 |
| Baseline LDH | 370 | 0.86 | 0.74, 0.96 | 0.019 | 0.87 | 0.75, 0.98 | 0.044 |
| ECOG 2+ | 370 | 0.51 | 0.31, 0.81 | 0.006 | 0.61 | 0.36, 1.03 | 0.068 |
| True EMD | 370 | 0.66 | 0.35, 1.20 | 0.2 | 0.62 | 0.31, 1.19 | 0.2 |
| Prior LOTs | 370 | 0.89 | 0.81, 0.97 | 0.014 | 0.95 | 0.85, 1.06 | 0.4 |
| Prior BCMA | 370 | 0.56 | 0.37, 0.86 | 0.009 | 0.52 | 0.31, 0.89 | 0.018 |
| Overall Survival | |||||||
| BsAb | 352 | ||||||
|     Tec | — | — | — | — | |||
| Â Â Â Â Elra | 1.41 | 0.96, 2.07 | 0.083 | 1.46 | 0.97, 2.20 | 0.073 | |
| Baseline Hgb | 352 | 0.73 | 0.66, 0.81 | <0.001 | 0.77 | 0.69, 0.86 | <0.001 |
| Baseline LDH | 352 | 1.13 | 1.06, 1.20 | <0.001 | 1.15 | 1.10, 1.20 | <0.001 |
| ECOG 2+ | 352 | 2.25 | 1.61, 3.16 | <0.001 | 1.92 | 1.32, 2.78 | <0.001 |
| True EMD | 352 | 1.34 | 0.88, 2.03 | 0.2 | 1.28 | 0.82, 2.00 | 0.3 |
| Prior LOTs | 352 | 1.06 | 0.99, 1.13 | 0.091 | 1.04 | 0.97, 1.11 | 0.3 |
| Prior BCMA | 352 | 1.09 | 0.77, 1.54 | 0.6 | 1.20 | 0.81, 1.78 | 0.4 |
| Progression-free Survival | |||||||
| BsAb | 355 | ||||||
|     Tec | — | — | — | — | |||
| Â Â Â Â Elra | 1.43 | 1.05, 1.94 | 0.021 | 1.52 | 1.11, 2.09 | 0.009 | |
| Baseline Hgb | 355 | 0.77 | 0.71, 0.83 | <0.001 | 0.81 | 0.74, 0.89 | <0.001 |
| Baseline LDH | 355 | 1.18 | 1.08, 1.28 | <0.001 | 1.17 | 1.08, 1.27 | <0.001 |
| ECOG 2+ | 355 | 1.63 | 1.23, 2.16 | <0.001 | 1.43 | 1.04, 1.95 | 0.027 |
| True EMD | 355 | 1.04 | 0.74, 1.47 | 0.8 | 0.94 | 0.67, 1.32 | 0.7 |
| Prior LOTs | 355 | 1.10 | 1.05, 1.16 | <0.001 | 1.06 | 1.00, 1.12 | 0.056 |
| Prior BCMA | 355 | 1.58 | 1.18, 2.10 | 0.002 | 1.63 | 1.17, 2.27 | 0.004 |
| Duration of Response | |||||||
| BsAb | 177 | ||||||
|     Tec | — | — | — | — | |||
| Â Â Â Â Elra | 1.49 | 0.87, 2.56 | 0.14 | 1.50 | 0.85, 2.63 | 0.2 | |
| Baseline Hgb | 177 | 0.84 | 0.74, 0.95 | 0.007 | 0.83 | 0.72, 0.96 | 0.014 |
| Baseline LDH | 177 | 1.05 | 0.90, 1.23 | 0.5 | 1.04 | 0.89, 1.21 | 0.7 |
| ECOG 2+ | 177 | 1.05 | 0.64, 1.72 | 0.8 | 0.89 | 0.52, 1.52 | 0.7 |
| True EMD | 177 | 1.48 | 0.75, 2.93 | 0.3 | 1.36 | 0.68, 2.72 | 0.4 |
| Prior LOTs | 177 | 1.01 | 0.89, 1.15 | 0.8 | 1.02 | 0.90, 1.17 | 0.7 |
| Prior BCMA | 177 | 1.08 | 0.67, 1.74 | 0.8 | 1.06 | 0.65, 1.75 | 0.8 |
| Infection-free Survival | |||||||
| BsAb | 358 | ||||||
|     Tec | — | — | — | — | |||
| Â Â Â Â Elra | 1.25 | 0.93, 1.67 | 0.13 | 1.24 | 0.92, 1.67 | 0.2 | |
| Baseline Hgb | 358 | 0.87 | 0.81, 0.94 | <0.001 | 0.90 | 0.83, 0.97 | 0.010 |
| Baseline LDH | 358 | 1.10 | 1.05, 1.15 | <0.001 | 1.09 | 1.05, 1.14 | <0.001 |
| ECOG 2+ | 358 | 1.70 | 1.30, 2.22 | <0.001 | 1.52 | 1.14, 2.02 | 0.004 |
| True EMD | 358 | 1.16 | 0.83, 1.62 | 0.4 | 1.10 | 0.79, 1.54 | 0.6 |
| Prior LOTs | 358 | 0.98 | 0.93, 1.03 | 0.4 | 0.98 | 0.93, 1.04 | 0.5 |
| Prior BCMA | 358 | 0.85 | 0.65, 1.11 | 0.2 | 0.99 | 0.73, 1.35 | >0.9 |
| ICANS 1+ | |||||||
| BsAb | 370 | ||||||
|     Tec | — | — | — | — | |||
| Â Â Â Â Elra | 0.87 | 0.49, 1.52 | 0.6 | 0.91 | 0.49, 1.65 | 0.8 | |
| Baseline Hgb | 370 | 0.83 | 0.72, 0.95 | 0.008 | 0.84 | 0.72, 0.98 | 0.024 |
| Baseline LDH | 370 | 1.05 | 0.95, 1.15 | 0.3 | 1.04 | 0.93, 1.14 | 0.5 |
| ECOG 2+ | 370 | 3.03 | 1.77, 5.21 | <0.001 | 2.44 | 1.39, 4.31 | 0.002 |
| True EMD | 370 | 1.40 | 0.69, 2.69 | 0.3 | 1.31 | 0.62, 2.62 | 0.5 |
| Prior LOTs | 370 | 0.88 | 0.77, 0.98 | 0.034 | 0.88 | 0.76, 1.01 | 0.069 |
| Prior BCMA | 370 | 0.61 | 0.36, 1.03 | 0.067 | 1.01 | 0.53, 1.92 | >0.9 |
| CRS 1+ | |||||||
| BsAb | 370 | ||||||
|     Tec | — | — | — | — | |||
| Â Â Â Â Elra | 0.53 | 0.34, 0.82 | 0.004 | 0.55 | 0.35, 0.86 | 0.010 | |
| Baseline Hgb | 370 | 0.94 | 0.85, 1.04 | 0.3 | 0.91 | 0.81, 1.01 | 0.076 |
| Baseline LDH | 370 | 0.98 | 0.90, 1.06 | 0.5 | 0.97 | 0.89, 1.06 | 0.5 |
| ECOG 2+ | 370 | 1.08 | 0.70, 1.68 | 0.7 | 0.97 | 0.60, 1.54 | 0.9 |
| True EMD | 370 | 0.58 | 0.33, 1.02 | 0.061 | 0.53 | 0.29, 0.94 | 0.033 |
| Prior LOTs | 370 | 0.92 | 0.84, 1.00 | 0.045 | 0.87 | 0.79, 0.96 | 0.006 |
| Prior BCMA | 370 | 1.17 | 0.77, 1.76 | 0.5 | 1.84 | 1.12, 3.04 | 0.016 |
| Steroid use | |||||||
| BsAb | 355 | ||||||
|     Tec | — | — | — | — | |||
| Â Â Â Â Elra | 0.83 | 0.48, 1.40 | 0.5 | 0.87 | 0.50, 1.50 | 0.6 | |
| Baseline Hgb | 355 | 0.89 | 0.79, 1.01 | 0.074 | 0.90 | 0.79, 1.03 | 0.14 |
| Baseline LDH | 355 | 1.07 | 0.98, 1.17 | 0.12 | 1.05 | 0.96, 1.15 | 0.3 |
| ECOG 2+ | 355 | 1.71 | 1.03, 2.83 | 0.036 | 1.32 | 0.78, 2.23 | 0.3 |
| True EMD | 355 | 1.61 | 0.85, 2.97 | 0.13 | 1.57 | 0.81, 2.94 | 0.2 |
| Prior LOTs | 355 | 0.88 | 0.79, 0.98 | 0.028 | 0.94 | 0.83, 1.06 | 0.3 |
| Prior BCMA | 355 | 0.47 | 0.29, 0.78 | 0.003 | 0.65 | 0.36, 1.16 | 0.15 |
| 1 Estimate = Odds ratio or hazard ratio | |||||||
data1 <- IPTW.data %>%
transmute(
site,
age_at_bs_ab_initiation,
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
bs_ab_name = factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra"), levels = c("Tec","Elra")),
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_bs_ab_step_up_d1,
date_of_pd,
date_of_last_contact,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
weights
)
dd <- datadist(data1)
options(datadist='dd')
ref_value = 85
S <- Surv(data1$Days.to.death.or.DLC, data1$Death)
f <- cph(S ~ rcs(age_at_bs_ab_initiation,3) * bs_ab_name , x=TRUE, y=TRUE,surv=TRUE,data=data1, weights = weights)
coxph( S ~ age_at_bs_ab_initiation * bs_ab_name , data=data1) %>% tbl_regression(exponentiate=TRUE)
| Characteristic | HR1 | 95% CI1 | p-value |
|---|---|---|---|
| age_at_bs_ab_initiation | 0.99 | 0.97, 1.01 | 0.3 |
| bs_ab_name | |||
|     Tec | — | — | |
| Â Â Â Â Elra | 0.05 | 0.00, 1.01 | 0.051 |
| age_at_bs_ab_initiation * bs_ab_name | |||
| Â Â Â Â age_at_bs_ab_initiation * Elra | 1.05 | 1.00, 1.09 | 0.029 |
| 1 HR = Hazard Ratio, CI = Confidence Interval | |||
model <- Predict(f, bs_ab_name, age_at_bs_ab_initiation, fun=function(x) exp(x) /exp(predict(f, data.frame(age_at_bs_ab_initiation = ref_value, bs_ab_name = "Tec")) ) )
OS_age <- ggplot(as.data.frame(model),aes(x=age_at_bs_ab_initiation, y=yhat, z = bs_ab_name)) +
geom_line(
aes(col = bs_ab_name)
) +
geom_ribbon(data = filter(model, bs_ab_name =="Elra"), aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0, fill = "steelblue" ) +
geom_ribbon(data = filter(model, bs_ab_name =="Tec"), aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0, fill = "coral2" ) +
theme_bw() +
labs(
title = "A)",
x = "Age (years)",
y = "OS (HR)") +
scale_color_manual(values = c ("coral2", "steelblue")) +
theme(title = element_text(size=16,face="bold"),legend.title=element_blank()) +
#coord_cartesian(ylim = c(0.5, 5.5)) +
geom_hline(yintercept = 1, linetype = "dashed", color = "gray", size = 0.5)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
OS_age
data1 <- IPTW.data %>%
transmute(
site,
age_at_bs_ab_initiation,
patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
bs_ab_name = factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra"), levels = c("Tec","Elra")),
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_bs_ab_step_up_d1,
date_of_pd,
date_of_last_contact,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1),
trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
weights
)
dd <- datadist(data1)
options(datadist='dd')
ref_value = 85
S <- Surv(data1$Days.to.death.or.relapse.or.DLC, data1$Relapse.or.death)
f <- cph(S ~ rcs(age_at_bs_ab_initiation,3) * bs_ab_name , x=TRUE, y=TRUE,surv=TRUE,data=data1, weights = weights)
coxph( S ~ age_at_bs_ab_initiation * bs_ab_name , data=data1) %>% tbl_regression(exponentiate=TRUE)
| Characteristic | HR1 | 95% CI1 | p-value |
|---|---|---|---|
| age_at_bs_ab_initiation | 0.98 | 0.97, 1.00 | 0.034 |
| bs_ab_name | |||
|     Tec | — | — | |
| Â Â Â Â Elra | 0.54 | 0.06, 4.75 | 0.6 |
| age_at_bs_ab_initiation * bs_ab_name | |||
| Â Â Â Â age_at_bs_ab_initiation * Elra | 1.01 | 0.98, 1.04 | 0.4 |
| 1 HR = Hazard Ratio, CI = Confidence Interval | |||
model <- Predict(f, bs_ab_name, age_at_bs_ab_initiation, fun=function(x) exp(x) /exp(predict(f, data.frame(age_at_bs_ab_initiation = ref_value, bs_ab_name = "Tec")) ) )
PFS_age <- ggplot(as.data.frame(model),aes(x=age_at_bs_ab_initiation, y=yhat, z = bs_ab_name)) +
geom_line(
aes(col = bs_ab_name)
) +
geom_ribbon(data = filter(model, bs_ab_name =="Elra"), aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0, fill = "steelblue" ) +
geom_ribbon(data = filter(model, bs_ab_name =="Tec"), aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0, fill = "coral2" ) +
theme_bw() +
labs(
title = "B)",
x = "Age (years)",
y = "PFS (HR)") +
scale_color_manual(values = c ("coral2", "steelblue")) +
theme(title = element_text(size=16,face="bold"),legend.title=element_blank()) +
#coord_cartesian(ylim = c(0.5, 5.5)) +
geom_hline(yintercept = 1, linetype = "dashed", color = "gray", size = 0.5)
PFS_age
combined_age <- (OS_age / PFS_age) +
plot_layout(guides = "collect", axes = "collect")
ggsave("svg/FIGURE5.svg", combined_age, device = "svg",
width = 6.5,
height = 6.5,
units = c("in"))
ggsave("jpeg/FIGURE5.jpeg", combined_age, device = "jpeg",
width = 6.5,
height = 6.5,
dpi = 600,
units = c("in"))
ggsave("tiff/FIGURE5.tiff", combined_age, device = "tiff",
width = 6.5,
height = 6.5,
dpi = 600,
units = c("in"),
compression = "lzw")
knitr::include_graphics("jpeg/FIGURE5.jpeg")
data1 <- IPTW.data %>%
filter(bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes==0) %>% ## Exclude patients who received bsAb as bridging
#filter(bs_ab_name == "tec") %>%
transmute(
bs_ab_name,
weights,
"CR/sCR" = ifelse(best_response %in% c("CR", "sCR"),1,0),
"PD" = ifelse(best_response %in% c("PD"),1,0),
"True EMD" = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2, 1,0),
"Age" = age_at_bs_ab_initiation/10,
"ECOG" = ecog_ps_at_bs_ab_initiation_0_5,
ECOG_2_4 = case_when(
is.na(ecog_ps_at_bs_ab_initiation_0_5) ~ NA,
ecog_ps_at_bs_ab_initiation_0_5 %in% c(2,3,4) ~ 1,
TRUE ~ 0),
"Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
"HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes | t_4_14_prior_to_bs_ab_0_no_1_yes | t_14_16_prior_to_bs_ab_0_no_1_yes | gain_amp1q21_prior_to_bs_ab_0_no_1_yes, 1, 0),
"Prior BCMA" = prior_bcma_directed_therapy_yes_no,
Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
Death = ifelse(death_yes_no,1,0),
date_of_pd,
date_of_last_contact,
date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death - date_of_bs_ab_step_up_d1),
"ALC BL" = baseline_alc,
"CrCl BL" = cr_cl_prior_to_step_up_dose_number_1,
"Trial ineligible" = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,0,1),
"Baseline platelets" = baseline_pl_ts,
baseline_hgb,
baseline_ldh = as.double(baseline_ldh_prior_to_bs_ab_initiation),
baseline_ferritin = as.numeric(baseline_ferritin_prior_to_bs_ab_initiation),
) %>% filter(!is.na(baseline_ldh))
## Tec
data_tec = data1 %>% filter(bs_ab_name == "tec")
dd_tec <- datadist(data_tec)
options(datadist='dd_tec')
S_tec <- Surv(data_tec$Days.to.death.or.relapse.or.DLC/30, data_tec$Relapse.or.death)
f_tec <- cph(S_tec ~ rcs(baseline_hgb, 3) * rcs(baseline_ldh, 3) , x=TRUE, y=TRUE,surv=TRUE,data=data_tec, weights=weights)
med_tec <- Quantile(f_tec)
pred_tec <- data.frame(Predict(f_tec, baseline_hgb, baseline_ldh,
fun = function(x) med_tec(lp = x))) %>%
mutate(drug = "Teclistamab")
## Elra
data_elra = data1 %>% filter(bs_ab_name == "elra")
dd_elra <- datadist(data_elra)
options(datadist='dd_elra')
S_elra <- Surv(data_elra$Days.to.death.or.relapse.or.DLC/30, data_elra$Relapse.or.death)
f_elra <- cph(S_elra ~ rcs(baseline_hgb, 3) * rcs(baseline_ldh, 3) , x=TRUE, y=TRUE,surv=TRUE,data=data_elra, weights=weights)
med_elra <- Quantile(f_elra)
pred_elra <- data.frame(Predict(f_elra, baseline_hgb, baseline_ldh,
fun = function(x) med_elra(lp = x))) %>%
mutate(drug = "Elranatamab")
## Combine both predictions
all_pred <- bind_rows(pred_tec, pred_elra)
z_range <- range(all_pred$yhat, na.rm = TRUE)
brks <- pretty(z_range, n = 10) # tweak n if you want more/fewer bands
all_pred <- all_pred %>%
mutate(
z_band = cut(
yhat,
breaks = brks,
include.lowest = TRUE,
right = FALSE
)
)
z_levels <- levels(all_pred$z_band)
cols <- viridis_pal()(length(z_levels))
contour_plot_tec <- all_pred %>%
filter(drug == "Teclistamab") %>%
ggplot(aes(baseline_hgb, baseline_ldh, fill = z_band)) +
geom_raster() + # or geom_tile(); blocky but guarantees alignment
theme_classic() +
theme(
panel.grid.major = element_line(color = "grey85"),
panel.grid.minor = element_line(color = "grey92")
)+
labs(
title = "A) Teclistamab",
x = "Baseline Hgb (g/dL)",
y = "Baseline LDH (U/L)",
fill = "Median PFS (months)"
) +
scale_fill_manual(
values = cols,
breaks = z_levels,
drop = FALSE,
na.translate = FALSE
) +
scale_x_continuous(breaks = seq(0, 20, 1))+
scale_y_continuous(breaks = seq(0, 2000, 100))+
coord_cartesian(xlim = c(7, 13), ylim = c(120, 780))
contour_plot_tec
contour_plot_elra <- all_pred %>%
filter(drug == "Elranatamab") %>%
ggplot(aes(baseline_hgb, baseline_ldh, fill = z_band)) +
geom_raster() +
theme_classic() +
theme(
legend.position = "none",
panel.grid.major = element_line(color = "grey85"),
panel.grid.minor = element_line(color = "grey92")
) +
labs(
title = "B) Elranatamab",
x = "Baseline Hgb (g/dL)",
y = "Baseline LDH (U/L)"
#fill = "Median PFS (months)"
) +
scale_fill_manual(
values = cols,
breaks = z_levels,
drop = FALSE,
na.translate = FALSE
) +
scale_x_continuous(breaks = seq(0, 20, 1))+
scale_y_continuous(breaks = seq(0, 2000, 100))+
coord_cartesian(xlim = c(7, 13), ylim = c(120, 780))
contour_plot_elra
combined_contour <- (contour_plot_tec / contour_plot_elra) +
plot_layout(guides = "collect", axes = "collect")
ggsave("svg/FIGURE6.svg", combined_contour, device = "svg",
width = 6.5,
height = 6.5,
units = c("in"))
ggsave("jpeg/FIGURE6.jpeg", combined_contour, device = "jpeg",
width = 6.5,
height = 6.5,
dpi = 600,
units = c("in"))
ggsave("tiff/FIGURE6.tiff", combined_contour, device = "tiff",
width = 6.5,
height = 6.5,
dpi = 600,
units = c("in"),
compression = "lzw")
knitr::include_graphics("jpeg/FIGURE6.jpeg")